import pandas as pd
import numpy as np
import os, glob, re, csv
import matplotlib.pyplot as plt
from scipy.signal import find_peaks, peak_widths, savgol_filter
from scipy.optimize import curve_fit
from matplotlib.tri import Triangulation
from matplotlib.colors import LinearSegmentedColormap, PowerNorm, Normalize
from pathlib import Path
from scipy.interpolate import griddata

# Define all deposition folders to process
BASE_DIR = r"C:\Dilan\Study\Master Thesis\Characterization\XRD"

DEPOSITIONS = [
    "30-09-2025 Dep 1",
    "22-09-2025 Dep 2",
    "24-09-2025 Dep 3",
    "07-10-2025 Dep 6",
    "08-10-2025 Dep 7",
    "09-10-2025 Dep 8"
]

# Target position (mm) - same for all depositions
target_x = 399.88
target_y = 224.1
target_z = 332.15

# Wafer center position (x=0, y=0 in XRD coordinates corresponds to this)
wafer_center_x = 242.5
wafer_center_y = 257.0
wafer_center_z = 435.0

print(f"Base directory: {BASE_DIR}")
print(f"\nDepositions to process: {len(DEPOSITIONS)}")
for dep in DEPOSITIONS:
    print(f"  - {dep}")

Base directory: C:\Dilan\Study\Master Thesis\Characterization\XRD

Depositions to process: 6
  - 30-09-2025 Dep 1
  - 22-09-2025 Dep 2
  - 24-09-2025 Dep 3
  - 07-10-2025 Dep 6
  - 08-10-2025 Dep 7
  - 09-10-2025 Dep 8

# Define XRD analysis parameters
REGIONS = [
    (20.5, 23.0),   # (100)
    (30.5, 33.0),   # (110)
    (37.5, 40.0),   # (111)
    (44.0, 47.0),   # (200)
    (50.0, 52.5),   # (210)
    (55.0, 57.5),   # (112)
]
REGION_NAMES = ["(100)", "(110)", "(111)", "(200)", "(210)", "(112)"]

# Peak finding settings
window_pts = 12
polyorder = 2
min_sep_deg = 0.25
width_deg = (0.06, 1.2)
ratio_min = 1.2
frac_min = 0.00

def _safe_window(n, w, p):
    if n < 5:
        return None, None
    w = min(w, n - (1 - (n % 2)))
    if w % 2 == 0: w += 1
    w = max(5, w)
    p = min(p, w - 2)
    return w, p

def extract_xy_from_filename(filename):
    """Extract X and Y coordinates from filename like 'something_X-value_Y-value.xy'"""
    name = os.path.basename(filename).rsplit('.', 1)[0]
    parts = name.split('_')
    try:
        x_str, y_str = parts[-2], parts[-1]
        def tok_to_float(t):
            if '-' in t[1:]:
                t = ('-' + t[1:].replace('-', '.', 1)) if t.startswith('-') else t.replace('-', '.', 1)
            return float(t)
        x_val, y_val = tok_to_float(x_str), tok_to_float(y_str)
        return x_val, y_val
    except Exception:
        return None, None

print("Analysis parameters configured.")

Analysis parameters configured.

# Process all depositions
all_deposition_results = {}

for dep_name in DEPOSITIONS:
    print(f"\n{'='*80}")
    print(f"Processing: {dep_name}")
    print(f"{'='*80}")
    
    dep_folder = os.path.join(BASE_DIR, dep_name)
    xy_data_folder = os.path.join(dep_folder, "XY Data")
    
    if not os.path.exists(xy_data_folder):
        print(f"  WARNING: XY Data folder not found in {dep_name}")
        continue
    
    os.chdir(xy_data_folder)
    files = sorted(glob.glob("*.xy"))
    
    if len(files) == 0:
        print(f"  WARNING: No .xy files found in {dep_name}")
        continue
    
    print(f"  Found {len(files)} XRD data files")
    
    # Initialize storage for this deposition
    all_file_results = []
    x_coordinates = []
    y_coordinates = []
    
    # Process each file
    for fp in files:
        data = np.loadtxt(fp)
        x, y = data[:,0], data[:,1]
        step = float(np.mean(np.diff(x)))
        
        # Smooth for peak detection
        w_use, p_use = _safe_window(len(x), window_pts, polyorder)
        ys = savgol_filter(y, w_use, p_use, mode="interp") if w_use is not None else y.copy()
        
        distance = max(1, int(min_sep_deg/step))
        width_s = (max(1, int(width_deg[0]/step)), max(2, int(width_deg[1]/step)))
        
        # Storage for THIS file: 6 values (one per region)
        file_fwhm = np.full(6, np.nan)
        file_intensity = np.full(6, np.nan)
        file_area = np.full(6, np.nan)
        file_peak_positions = [None] * 6
        
        per_region_count = {i: 0 for i in range(1, len(REGIONS)+1)}
        
        for rid, (lo, hi) in enumerate(REGIONS, start=1):
            m = (x >= lo) & (x <= hi)
            if np.count_nonzero(m) < 7:
                continue
            
            xr = x[m]
            yr = ys[m]
            idx0 = np.nonzero(m)[0][0]
            
            # Find peaks
            cand_r, prop = find_peaks(yr, distance=distance, width=width_s)
            if cand_r.size == 0:
                continue
            
            keep_r = []
            for i, pk in enumerate(cand_r):
                lb = prop["left_bases"][i]
                rb = prop["right_bases"][i]
                base = max(yr[lb], yr[rb])
                amp = yr[pk] - base
                if amp <= 0: continue
                if yr[pk] / max(base, 1e-9) < ratio_min: continue
                if amp < frac_min * base: continue
                keep_r.append(i)
            
            if not keep_r:
                continue
            
            peaks_r = cand_r[keep_r]
            widths, w_left, w_right, _ = peak_widths(yr, peaks_r, rel_height=0.5)
            
            idx_axis = np.arange(len(xr), dtype=float)
            left_x = np.interp(w_left, idx_axis, xr)
            right_x = np.interp(w_right, idx_axis, xr)
            fwhm_deg = np.abs(right_x - left_x)
            
            for j, pk in enumerate(peaks_r):
                lb = prop["left_bases"][keep_r[j]]
                rb = prop["right_bases"][keep_r[j]]
                
                x_lb, y_lb = xr[lb], yr[lb]
                x_rb, y_rb = xr[rb], yr[rb]
                
                lx, rx = x_lb, x_rb
                if lx >= rx:
                    lx = left_x[j]
                    rx = right_x[j]
                
                seg_mask = (xr >= lx) & (xr <= rx)
                if np.count_nonzero(seg_mask) < 2:
                    k0 = max(0, pk - 3)
                    k1 = min(len(xr), pk + 4)
                    seg_mask = np.zeros_like(xr, dtype=bool)
                    seg_mask[k0:k1] = True
                
                xs = xr[seg_mask]
                ys_seg = yr[seg_mask]
                
                if x_rb != x_lb:
                    slope = (y_rb - y_lb) / (x_rb - x_lb)
                    y_base = y_lb + slope * (xs - x_lb)
                else:
                    y_base = np.full_like(xs, np.mean([y_lb, y_rb]))
                
                area = np.trapz(np.maximum(ys_seg - y_base, 0.0), xs)
                
                if area == 0.0 or not np.isfinite(area):
                    peak_height = yr[pk] - np.mean([y_lb, y_rb])
                    area = peak_height * fwhm_deg[j] * 0.5 * np.sqrt(np.pi / np.log(2))
                
                gi = idx0 + pk
                per_region_count[rid] += 1
                
                # Store first peak per region
                if per_region_count[rid] == 1:
                    file_fwhm[rid-1] = float(fwhm_deg[j])
                    file_intensity[rid-1] = float(yr[pk])
                    file_area[rid-1] = float(area)
                    file_peak_positions[rid-1] = float(x[gi])
        
        # Extract X and Y coordinates from filename
        x_val, y_val = extract_xy_from_filename(fp)
        if x_val is not None and y_val is not None:
            x_coordinates.append(x_val)
            y_coordinates.append(y_val)
        
        # Store results for this file
        all_file_results.append({
            'fwhm': file_fwhm,
            'intensity': file_intensity,
            'area': file_area,
            'peak_positions': file_peak_positions,
            'filename': fp
        })
    
    # Convert to numpy arrays
    results = np.array([f['fwhm'] for f in all_file_results])
    peak_intensity = np.array([f['intensity'] for f in all_file_results])
    peak_areas = np.array([f['area'] for f in all_file_results])
    
    # Take absolute values of FWHM
    results = np.abs(results)
    
    # Store results for this deposition
    all_deposition_results[dep_name] = {
        'x_coordinates': np.array(x_coordinates),
        'y_coordinates': np.array(y_coordinates),
        'fwhm': results,
        'peak_intensity': peak_intensity,
        'peak_areas': peak_areas,
        'num_files': len(files),
        'num_processed': len(all_file_results),
        'files': files,
        'file_results': all_file_results
    }
    
    print(f"  Processed {len(all_file_results)} files successfully")

print(f"\n\n{'='*80}")
print(f"Processing complete! Results stored for {len(all_deposition_results)} depositions.")
print(f"{'='*80}")

================================================================================
Processing: 30-09-2025 Dep 1
================================================================================
  Found 61 XRD data files
  Processed 61 files successfully

================================================================================
Processing: 22-09-2025 Dep 2
================================================================================
  Found 71 XRD data files
  Processed 71 files successfully

================================================================================
Processing: 24-09-2025 Dep 3
================================================================================
  Found 61 XRD data files
  Processed 61 files successfully

================================================================================
Processing: 07-10-2025 Dep 6
================================================================================
  Found 61 XRD data files
  Processed 61 files successfully

================================================================================
Processing: 08-10-2025 Dep 7
================================================================================
  Found 61 XRD data files
  Processed 61 files successfully

================================================================================
Processing: 09-10-2025 Dep 8
================================================================================
  Found 61 XRD data files
  Processed 61 files successfully


================================================================================
Processing complete! Results stored for 6 depositions.
================================================================================

# Create stacked plots for each deposition
for dep_name in DEPOSITIONS:
    if dep_name not in all_deposition_results:
        continue
    
    print(f"\nCreating stacked plot for {dep_name}")
    
    dep_folder = os.path.join(BASE_DIR, dep_name)
    xy_data_folder = os.path.join(dep_folder, "XY Data")
    os.chdir(xy_data_folder)
    
    files = all_deposition_results[dep_name]['files']
    
    # Select files to plot (first, middle points, last - up to 9 files)
    if len(files) <= 9:
        files_to_plot = files
    else:
        # Select evenly spaced files
        indices = np.linspace(0, len(files)-1, 9, dtype=int)
        files_to_plot = [files[i] for i in indices]
    
    # Create ONE figure for all plots
    fig, ax = plt.subplots(figsize=(12, 14))
    
    # Define a colormap
    colors = plt.cm.tab10(np.linspace(0, 1, len(files_to_plot)))
    
    # Load data and prepare labels
    data_list = []
    label_list = []
    
    for f in files_to_plot:
        df = pd.read_csv(f, delim_whitespace=True, header=None, names=["X", "Y"])
        data_list.append(df)
        
        # Extract coordinates and create label
        x_coord, y_coord = extract_xy_from_filename(f)
        if x_coord is not None and y_coord is not None:
            label = f"BTO Mapped X: {x_coord:.1f}, Y: {y_coord:.1f}"
        else:
            label = os.path.basename(f)
        label_list.append(label)
    
    # Use logarithmic offset: multiply each curve by 10^idx
    log_offset_factor = 10
    
    # Store handles and labels for reverse ordering
    handles = []
    labels_for_legend = []
    
    # Plot with logarithmic offsets
    for idx, (df, label) in enumerate(zip(data_list, label_list)):
        offset_multiplier = log_offset_factor ** idx
        
        line, = ax.plot(df["X"], df["Y"] * offset_multiplier,
                        color=colors[idx], linewidth=1.2)
        
        handles.append(line)
        labels_for_legend.append(label)
        
        # Optional: Add horizontal gridline at baseline of each curve
        baseline = df["Y"].min() * offset_multiplier
        ax.axhline(y=baseline, color=colors[idx], linestyle='--',
                   alpha=0.2, linewidth=0.5)
    
    # Set logarithmic scale for Y-axis
    ax.set_yscale('log')
    
    # Reverse the order for legend
    handles_reversed = handles[::-1]
    labels_reversed = labels_for_legend[::-1]
    
    # Add legend outside the plot on the right
    ax.legend(handles_reversed, labels_reversed, 
              loc='center left', fontsize=11, framealpha=0.95,
              bbox_to_anchor=(1.02, 0.5))
    
    # Add labels and formatting
    ax.set_title(f"{dep_name}: Stacked XRD Data from Different Locations (Log Scale)", fontsize=14)
    ax.set_xlabel("2θ (deg)", fontsize=12)
    ax.set_ylabel("Intensity (a.u.) - log scale with offset", fontsize=12)
    ax.set_xlim(20, 70)
    ax.grid(True, alpha=0.3, axis='x')
    ax.grid(True, alpha=0.2, axis='y', which='both')
    plt.tight_layout()
    plt.show()

Creating stacked plot for 30-09-2025 Dep 1

Creating stacked plot for 22-09-2025 Dep 2

Creating stacked plot for 24-09-2025 Dep 3

Creating stacked plot for 07-10-2025 Dep 6

Creating stacked plot for 08-10-2025 Dep 7

Creating stacked plot for 09-10-2025 Dep 8

# Create individual plots with peak detection for each deposition
for dep_name in DEPOSITIONS:
    if dep_name not in all_deposition_results:
        continue
    
    print(f"\n{'='*80}")
    print(f"Individual plots for: {dep_name}")
    print(f"{'='*80}")
    
    dep_folder = os.path.join(BASE_DIR, dep_name)
    xy_data_folder = os.path.join(dep_folder, "XY Data")
    os.chdir(xy_data_folder)
    
    files = all_deposition_results[dep_name]['files']
    file_results = all_deposition_results[dep_name]['file_results']
    
    # Process each file and create individual plot
    for fp, file_res in zip(files, file_results):
        data = np.loadtxt(fp)
        x, y = data[:,0], data[:,1]
        step = float(np.mean(np.diff(x)))
        
        # Smooth
        w_use, p_use = _safe_window(len(x), window_pts, polyorder)
        ys = savgol_filter(y, w_use, p_use, mode="interp") if w_use is not None else y.copy()
        
        distance = max(1, int(min_sep_deg/step))
        width_s = (max(1, int(width_deg[0]/step)), max(2, int(width_deg[1]/step)))
        
        all_peaks_idx = []
        all_fwhm_deg = []
        all_area = []
        all_region_id = []
        all_Ipk = []
        all_labels = []
        
        per_region_count = {i: 0 for i in range(1, len(REGIONS)+1)}
        
        for rid, (lo, hi) in enumerate(REGIONS, start=1):
            m = (x >= lo) & (x <= hi)
            if np.count_nonzero(m) < 7:
                continue
            
            xr = x[m]
            yr = ys[m]
            idx0 = np.nonzero(m)[0][0]
            
            cand_r, prop = find_peaks(yr, distance=distance, width=width_s)
            if cand_r.size == 0:
                continue
            
            keep_r = []
            for i, pk in enumerate(cand_r):
                lb = prop["left_bases"][i]
                rb = prop["right_bases"][i]
                base = max(yr[lb], yr[rb])
                amp = yr[pk] - base
                if amp <= 0: continue
                if yr[pk] / max(base, 1e-9) < ratio_min: continue
                if amp < frac_min * base: continue
                keep_r.append(i)
            
            if not keep_r:
                continue
            
            peaks_r = cand_r[keep_r]
            widths, w_left, w_right, _ = peak_widths(yr, peaks_r, rel_height=0.5)
            
            idx_axis = np.arange(len(xr), dtype=float)
            left_x = np.interp(w_left, idx_axis, xr)
            right_x = np.interp(w_right, idx_axis, xr)
            fwhm_deg = np.abs(right_x - left_x)
            
            for j, pk in enumerate(peaks_r):
                lb = prop["left_bases"][keep_r[j]]
                rb = prop["right_bases"][keep_r[j]]
                
                x_lb, y_lb = xr[lb], yr[lb]
                x_rb, y_rb = xr[rb], yr[rb]
                
                lx, rx = x_lb, x_rb
                if lx >= rx:
                    lx = left_x[j]
                    rx = right_x[j]
                
                seg_mask = (xr >= lx) & (xr <= rx)
                if np.count_nonzero(seg_mask) < 2:
                    k0 = max(0, pk - 3)
                    k1 = min(len(xr), pk + 4)
                    seg_mask = np.zeros_like(xr, dtype=bool)
                    seg_mask[k0:k1] = True
                
                xs = xr[seg_mask]
                ys_seg = yr[seg_mask]
                
                if x_rb != x_lb:
                    slope = (y_rb - y_lb) / (x_rb - x_lb)
                    y_base = y_lb + slope * (xs - x_lb)
                else:
                    y_base = np.full_like(xs, np.mean([y_lb, y_rb]))
                
                area = np.trapz(np.maximum(ys_seg - y_base, 0.0), xs)
                
                if area == 0.0 or not np.isfinite(area):
                    peak_height = yr[pk] - np.mean([y_lb, y_rb])
                    area = peak_height * fwhm_deg[j] * 0.5 * np.sqrt(np.pi / np.log(2))
                
                gi = idx0 + pk
                all_peaks_idx.append(gi)
                all_fwhm_deg.append(float(fwhm_deg[j]))
                all_area.append(float(area))
                all_region_id.append(rid)
                all_Ipk.append(float(yr[pk]))
                
                # Name the peak
                base_name = REGION_NAMES[rid-1] if rid-1 < len(REGION_NAMES) else f"R{rid}"
                per_region_count[rid] += 1
                label = base_name if per_region_count[rid] == 1 else f"{base_name}-{per_region_count[rid]}"
                all_labels.append(label)
        
        # Print peak info
        print(f"\n{os.path.basename(fp)}  peaks: {len(all_peaks_idx)}")
        for lbl, rid, gi, Ipk, fwhm, area in zip(all_labels, all_region_id, all_peaks_idx, all_Ipk, all_fwhm_deg, all_area):
            print(f"  {lbl:8s} [R{rid}]  {x[gi]:7.3f}°   I≈{Ipk:.1f}   FWHM≈{fwhm:.4f}°   Area≈{area:.2f}")
        
        # Plot
        plt.figure(figsize=(7.6, 4.8))
        plt.plot(x, ys, label="smoothed", color="0.2")
        
        # Annotate peaks
        if all_peaks_idx:
            p_x = x[np.array(all_peaks_idx)]
            p_y = np.array(all_Ipk)
            plt.plot(p_x, p_y, "rx", label="peaks")
            for (xx, yy, lbl) in zip(p_x, p_y, all_labels):
                plt.annotate(lbl, (xx, yy), xytext=(0, 8), textcoords="offset points",
                             ha="center", va="bottom", fontsize=8, rotation=0,
                             bbox=dict(boxstyle="round,pad=0.2", fc="white", ec="none", alpha=0.6))
        
        # Draw region spans
        for (lo, hi), nm in zip(REGIONS, REGION_NAMES):
            plt.axvspan(lo, hi, alpha=0.08)
            mid = 0.5*(lo+hi)
            plt.text(mid, plt.ylim()[1]*0.85, nm, ha="center", va="top", fontsize=8, alpha=0.6)
        
        # Extract coordinates for title
        x_val, y_val = extract_xy_from_filename(fp)
        if x_val is not None and y_val is not None:
            title_text = f"{dep_name} — X: {x_val:.2f}  Y: {y_val:.2f}"
        else:
            title_text = f"{dep_name} — {os.path.basename(fp)}"
        
        plt.title(title_text)
        plt.xlabel("2θ (deg)")
        plt.ylabel("Intensity (a.u.)")
        plt.yscale("log")
        plt.legend()
        plt.tight_layout()
        plt.show()

================================================================================
Individual plots for: 30-09-2025 Dep 1
================================================================================

BTO_Mapping_Theta2Theta_30-09-2025_001_01_0-000_0-000.xy  peaks: 6
  (100)    [R1]   22.050°   I≈814.6   FWHM≈1.1210°   Area≈136.91
  (110)    [R2]   31.370°   I≈3680.8   FWHM≈1.7841°   Area≈1169.32
  (111)    [R3]   38.730°   I≈843.7   FWHM≈1.4478°   Area≈224.16
  (200)    [R4]   45.070°   I≈462.6   FWHM≈2.1301°   Area≈139.24
  (210)    [R5]   50.770°   I≈280.7   FWHM≈1.7390°   Area≈57.95
  (112)    [R6]   55.960°   I≈658.6   FWHM≈1.7478°   Area≈231.85

BTO_Mapping_Theta2Theta_30-09-2025_001_02_10-000_0-000.xy  peaks: 6
  (100)    [R1]   22.040°   I≈864.5   FWHM≈1.1144°   Area≈122.19
  (110)    [R2]   31.360°   I≈3844.4   FWHM≈1.7868°   Area≈1210.98
  (111)    [R3]   38.710°   I≈812.8   FWHM≈1.4590°   Area≈214.00
  (200)    [R4]   45.060°   I≈496.2   FWHM≈2.1197°   Area≈147.90
  (210)    [R5]   50.720°   I≈295.7   FWHM≈1.8346°   Area≈63.45
  (112)    [R6]   55.920°   I≈674.2   FWHM≈1.7657°   Area≈248.42

BTO_Mapping_Theta2Theta_30-09-2025_001_03_20-000_0-000.xy  peaks: 6
  (100)    [R1]   22.060°   I≈831.9   FWHM≈1.1378°   Area≈131.63
  (110)    [R2]   31.370°   I≈4053.6   FWHM≈1.7873°   Area≈1269.06
  (111)    [R3]   38.710°   I≈801.5   FWHM≈1.4581°   Area≈233.07
  (200)    [R4]   45.050°   I≈487.2   FWHM≈2.1249°   Area≈146.10
  (210)    [R5]   50.710°   I≈292.0   FWHM≈1.8648°   Area≈62.64
  (112)    [R6]   55.940°   I≈670.2   FWHM≈1.7728°   Area≈260.74

BTO_Mapping_Theta2Theta_30-09-2025_001_04_30-000_0-000.xy  peaks: 6
  (100)    [R1]   22.040°   I≈830.9   FWHM≈1.1376°   Area≈166.32
  (110)    [R2]   31.360°   I≈4132.1   FWHM≈1.7896°   Area≈1319.48
  (111)    [R3]   38.720°   I≈813.1   FWHM≈1.4541°   Area≈227.66
  (200)    [R4]   45.050°   I≈476.7   FWHM≈2.1278°   Area≈149.51
  (210)    [R5]   50.660°   I≈294.8   FWHM≈1.8290°   Area≈67.21
  (112)    [R6]   55.910°   I≈675.9   FWHM≈1.7746°   Area≈254.63

BTO_Mapping_Theta2Theta_30-09-2025_001_05_40-000_0-000.xy  peaks: 6
  (100)    [R1]   22.070°   I≈805.4   FWHM≈1.1309°   Area≈119.36
  (110)    [R2]   31.370°   I≈4233.9   FWHM≈1.7878°   Area≈1349.20
  (111)    [R3]   38.730°   I≈799.0   FWHM≈1.4555°   Area≈237.61
  (200)    [R4]   45.080°   I≈444.4   FWHM≈2.1294°   Area≈127.37
  (210)    [R5]   50.660°   I≈293.5   FWHM≈1.8664°   Area≈61.65
  (112)    [R6]   55.960°   I≈674.4   FWHM≈1.7698°   Area≈262.26

BTO_Mapping_Theta2Theta_30-09-2025_001_06_9-135_4-067.xy  peaks: 6
  (100)    [R1]   22.040°   I≈816.6   FWHM≈1.1316°   Area≈112.28
  (110)    [R2]   31.360°   I≈3728.8   FWHM≈1.7787°   Area≈1161.11
  (111)    [R3]   38.730°   I≈861.4   FWHM≈1.4439°   Area≈231.78
  (200)    [R4]   45.060°   I≈449.3   FWHM≈2.0967°   Area≈123.74
  (210)    [R5]   50.690°   I≈278.1   FWHM≈1.7331°   Area≈59.12
  (112)    [R6]   55.940°   I≈653.9   FWHM≈1.7580°   Area≈240.85

BTO_Mapping_Theta2Theta_30-09-2025_001_07_18-271_8-135.xy  peaks: 6
  (100)    [R1]   22.090°   I≈760.7   FWHM≈1.1330°   Area≈100.13
  (110)    [R2]   31.380°   I≈3689.2   FWHM≈1.7799°   Area≈1172.76
  (111)    [R3]   38.710°   I≈905.6   FWHM≈1.4381°   Area≈256.23
  (200)    [R4]   45.070°   I≈417.1   FWHM≈2.1072°   Area≈106.95
  (210)    [R5]   50.700°   I≈280.7   FWHM≈1.7323°   Area≈63.11
  (112)    [R6]   55.920°   I≈664.3   FWHM≈1.7456°   Area≈247.53

BTO_Mapping_Theta2Theta_30-09-2025_001_08_27-406_12-202.xy  peaks: 6
  (100)    [R1]   22.050°   I≈711.6   FWHM≈1.1318°   Area≈88.79
  (110)    [R2]   31.390°   I≈3618.5   FWHM≈1.7750°   Area≈1141.72
  (111)    [R3]   38.740°   I≈952.0   FWHM≈1.4319°   Area≈265.79
  (200)    [R4]   45.130°   I≈375.4   FWHM≈2.0642°   Area≈110.47
  (210)    [R5]   50.720°   I≈276.6   FWHM≈1.7098°   Area≈53.30
  (112)    [R6]   55.960°   I≈642.9   FWHM≈1.7456°   Area≈246.35

BTO_Mapping_Theta2Theta_30-09-2025_001_09_36-542_16-269.xy  peaks: 6
  (100)    [R1]   22.100°   I≈696.9   FWHM≈1.0866°   Area≈81.56
  (110)    [R2]   31.390°   I≈3507.8   FWHM≈1.7665°   Area≈1117.34
  (111)    [R3]   38.720°   I≈925.7   FWHM≈1.4222°   Area≈265.51
  (200)    [R4]   45.130°   I≈330.9   FWHM≈1.8226°   Area≈80.84
  (210)    [R5]   50.670°   I≈266.9   FWHM≈1.6759°   Area≈55.22
  (112)    [R6]   56.010°   I≈639.4   FWHM≈1.7355°   Area≈249.24

BTO_Mapping_Theta2Theta_30-09-2025_001_10_6-691_7-431.xy  peaks: 6
  (100)    [R1]   22.050°   I≈771.5   FWHM≈1.1032°   Area≈102.88
  (110)    [R2]   31.380°   I≈3605.3   FWHM≈1.7734°   Area≈1118.03
  (111)    [R3]   38.710°   I≈902.3   FWHM≈1.4380°   Area≈257.85
  (200)    [R4]   45.090°   I≈414.1   FWHM≈2.0937°   Area≈119.17
  (210)    [R5]   50.680°   I≈280.0   FWHM≈1.7632°   Area≈58.59
  (112)    [R6]   55.940°   I≈658.3   FWHM≈1.7405°   Area≈240.05

# Helper functions for contour plotting
def _unique_xy(x, y, z):
    arr = np.column_stack([x, y, z])
    keys = np.ascontiguousarray(arr[:, :2]).view(
        np.dtype((np.void, arr[:, :2].dtype.itemsize * 2))
    ).ravel()
    _, idx = np.unique(keys, return_index=True)
    idx = np.sort(idx)
    return arr[idx, 0], arr[idx, 1], arr[idx, 2]

def _contour_one(ax, xs, ys, zs, title, nlevels=100, cmap="coolwarm", try_grid_fallback=True):
    m = np.isfinite(xs) & np.isfinite(ys) & np.isfinite(zs)
    xs, ys, zs = xs[m], ys[m], zs[m]
    if len(xs) == 0:
        ax.text(0.5, 0.5, "No data", ha='center', va='center', transform=ax.transAxes)
        ax.set_axis_off()
        return None

    xs, ys, zs = _unique_xy(xs, ys, zs)
    if len(xs) < 3:
        ax.text(0.5, 0.5, "Too few points", ha='center', va='center', transform=ax.transAxes)
        ax.set_axis_off()
        return None

    # Robust color limits
    vmin = np.nanpercentile(zs, 5)
    vmax = np.nanpercentile(zs, 95)
    if not np.isfinite(vmin) or not np.isfinite(vmax) or vmin == vmax:
        vmin, vmax = float(np.nanmin(zs)), float(np.nanmax(zs))
        if not np.isfinite(vmin) or vmin == vmax:
            ax.text(0.5, 0.5, "No variation", ha='center', va='center', transform=ax.transAxes)
            ax.set_axis_off()
            return None

    # Triangulation or grid fallback
    try:
        tri = Triangulation(xs, ys)
        cf = ax.tricontourf(tri, zs, levels=np.linspace(vmin, vmax, nlevels),
                            extend='both', cmap=cmap)
    except Exception:
        if not try_grid_fallback:
            ax.text(0.5, 0.5, "Triangulation failed", ha='center', va='center', transform=ax.transAxes)
            ax.set_axis_off()
            return None

        gx = np.linspace(xs.min(), xs.max(), 60)
        gy = np.linspace(ys.min(), ys.max(), 60)
        GX, GY = np.meshgrid(gx, gy)
        GZ = griddata((xs, ys), zs, (GX, GY), method='linear')
        if np.isnan(GZ).all():
            GZ = griddata((xs, ys), zs, (GX, GY), method='nearest')
        if np.isnan(GZ).all():
            ax.text(0.5, 0.5, "Insufficient geometry", ha='center', va='center', transform=ax.transAxes)
            ax.set_axis_off()
            return None

        cf = ax.contourf(GX, GY, GZ, levels=np.linspace(vmin, vmax, nlevels),
                         extend='both', cmap=cmap)

    ax.scatter(xs, ys, s=10, c='k', alpha=0.3)
    ax.set_title(title, fontsize=11)
    ax.set_xlabel("X [mm]")
    ax.set_ylabel("Y [mm]")
    return cf

# Normalize data into shape (N, 6)
def normalize_six(arr, N):
    a = np.asarray(arr)
    if a.ndim == 2 and a.shape == (N, 6):
        return a
    if a.ndim == 1 and a.size == N * 6:
        return a.reshape(N, 6)
    if a.ndim == 2 and a.shape[0] == N * 6 and a.shape[1] == 6:
        blocks = np.split(a, 6, axis=0)
        return np.column_stack([blocks[j][:, j] for j in range(6)])
    if a.ndim == 2 and a.shape == (N * 6, 1):
        return a.reshape(N, 6)
    raise ValueError(f"Cannot reshape array with shape {a.shape} to (N,6).")

# Create FWHM contour maps for all depositions
for dep_name in DEPOSITIONS:
    if dep_name not in all_deposition_results:
        continue
    
    print(f"\nCreating FWHM contour maps for {dep_name}")
    
    x_pos = all_deposition_results[dep_name]['x_coordinates']
    y_pos = all_deposition_results[dep_name]['y_coordinates']
    fwhm = all_deposition_results[dep_name]['fwhm']
    N = x_pos.size
    
    fwhm = normalize_six(fwhm, N)
    
    titles = [f"Region {REGION_NAMES[i]}" for i in range(6)]
    fig, axes = plt.subplots(2, 3, figsize=(14, 8), constrained_layout=True)
    cfs = []
    for i in range(6):
        r, c = divmod(i, 3)
        cf = _contour_one(axes[r, c], x_pos, y_pos, fwhm[:, i],
                          title=f"{titles[i]}: FWHM (°)", cmap="coolwarm")
        cfs.append(cf)
    
    # Row-wise colorbars
    for r in range(2):
        row = [cfs[r*3 + j] for j in range(3) if cfs[r*3 + j] is not None]
        if row:
            fig.colorbar(row[0], ax=axes[r, :], shrink=0.95, label="FWHM (°)")
    
    plt.suptitle(f"{dep_name}: Contour maps of FWHM per region", y=1.02, fontsize=14)
    plt.show()

Creating FWHM contour maps for 30-09-2025 Dep 1

Creating FWHM contour maps for 22-09-2025 Dep 2

Creating FWHM contour maps for 24-09-2025 Dep 3

Creating FWHM contour maps for 07-10-2025 Dep 6

Creating FWHM contour maps for 08-10-2025 Dep 7

Creating FWHM contour maps for 09-10-2025 Dep 8

# Create Peak Intensity contour maps for all depositions
for dep_name in DEPOSITIONS:
    if dep_name not in all_deposition_results:
        continue
    
    print(f"\nCreating Peak Intensity contour maps for {dep_name}")
    
    x_pos = all_deposition_results[dep_name]['x_coordinates']
    y_pos = all_deposition_results[dep_name]['y_coordinates']
    peakI = all_deposition_results[dep_name]['peak_intensity']
    N = x_pos.size
    
    peakI = normalize_six(peakI, N)
    
    titles = [f"Region {REGION_NAMES[i]}" for i in range(6)]
    fig, axes = plt.subplots(2, 3, figsize=(14, 8), constrained_layout=True)
    cfs = []
    for i in range(6):
        r, c = divmod(i, 3)
        cf = _contour_one(axes[r, c], x_pos, y_pos, peakI[:, i],
                          title=f"{titles[i]}: Peak intensity [counts per second]",
                          cmap="coolwarm", nlevels=120)
        cfs.append(cf)
    
    # Row-wise colorbars
    for r in range(2):
        row = [cfs[r*3 + j] for j in range(3) if cfs[r*3 + j] is not None]
        if row:
            fig.colorbar(row[0], ax=axes[r, :], shrink=0.95, label="Peak intensity")
    
    plt.suptitle(f"{dep_name}: Contour maps of peak intensity per region", y=1.02, fontsize=14)
    plt.show()

Creating Peak Intensity contour maps for 30-09-2025 Dep 1

Creating Peak Intensity contour maps for 22-09-2025 Dep 2

Creating Peak Intensity contour maps for 24-09-2025 Dep 3

Creating Peak Intensity contour maps for 07-10-2025 Dep 6

Creating Peak Intensity contour maps for 08-10-2025 Dep 7

Creating Peak Intensity contour maps for 09-10-2025 Dep 8

# Helper function for contour with custom levels
def _robust_minmax(zs):
    vmin = np.nanpercentile(zs, 5)
    vmax = np.nanpercentile(zs, 95)
    if not np.isfinite(vmin) or not np.isfinite(vmax) or vmin == vmax:
        vmin = float(np.nanmin(zs))
        vmax = float(np.nanmax(zs))
    return vmin, vmax

def _contour_one_custom(ax, xs_all, ys_all, zs_all, title, nlevels=100, cmap="coolwarm", levels=None):
    """Create smooth contour plot using triangulation. Shows all measurement points with markers."""
    mask_finite = np.isfinite(xs_all) & np.isfinite(ys_all) & np.isfinite(zs_all)
    xs = xs_all[mask_finite]
    ys = ys_all[mask_finite]
    zs = zs_all[mask_finite]
    
    if len(xs) == 0:
        ax.text(0.5, 0.5, "No data", ha="center", va="center", transform=ax.transAxes)
        ax.scatter(xs_all, ys_all, s=30, c='red', marker='x', alpha=0.6, linewidths=2)
        ax.set_title(title, fontsize=11)
        ax.set_xlabel("X [mm]")
        ax.set_ylabel("Y [mm]")
        ax.set_aspect('equal', adjustable='box')
        return None

    xs, ys, zs = _unique_xy(xs, ys, zs)
    
    if len(xs) < 3:
        ax.text(0.5, 0.5, "Too few points", ha="center", va="center", transform=ax.transAxes)
        ax.scatter(xs_all, ys_all, s=30, c='red', marker='x', alpha=0.6, linewidths=2)
        ax.set_title(title, fontsize=11)
        ax.set_xlabel("X [mm]")
        ax.set_ylabel("Y [mm]")
        ax.set_aspect('equal', adjustable='box')
        return None

    finite_levels = None
    if levels is not None:
        finite_levels = np.asarray(levels, dtype=float)
        finite_levels = finite_levels[np.isfinite(finite_levels)]

    if finite_levels is None or finite_levels.size == 0:
        vmin, vmax = _robust_minmax(zs)
        if not np.isfinite(vmin) or not np.isfinite(vmax) or vmin == vmax:
            ax.text(0.5, 0.5, "No variation", ha="center", va="center", transform=ax.transAxes)
            ax.set_title(title, fontsize=11)
            ax.set_xlabel("X [mm]")
            ax.set_ylabel("Y [mm]")
            ax.set_aspect('equal', adjustable='box')
            return None
        finite_levels = np.linspace(vmin, vmax, nlevels)

    try:
        tri = Triangulation(xs, ys)
        cf = ax.tricontourf(tri, zs, levels=finite_levels, extend="both", cmap=cmap)
    except Exception:
        gx = np.linspace(xs.min(), xs.max(), 100)
        gy = np.linspace(ys.min(), ys.max(), 100)
        GX, GY = np.meshgrid(gx, gy)
        GZ = griddata((xs, ys), zs, (GX, GY), method="linear")
        if np.isnan(GZ).all():
            GZ = griddata((xs, ys), zs, (GX, GY), method="nearest")
        if np.isnan(GZ).all():
            ax.text(0.5, 0.5, "Triangulation failed", ha="center", va="center", transform=ax.transAxes)
            ax.set_title(title, fontsize=11)
            ax.set_xlabel("X [mm]")
            ax.set_ylabel("Y [mm]")
            ax.set_aspect('equal', adjustable='box')
            return None
        
        cf = ax.contourf(GX, GY, GZ, levels=finite_levels, extend="both", cmap=cmap)

    # Scatter plot: show ALL measurement points
    if np.sum(mask_finite) > 0:
        ax.scatter(xs_all[mask_finite], ys_all[mask_finite], 
                   s=30, c="k", alpha=0.6, edgecolors='white', linewidths=0.5, zorder=5)
    
    if np.sum(~mask_finite) > 0:
        ax.scatter(xs_all[~mask_finite], ys_all[~mask_finite], 
                   s=30, c="red", marker='x', alpha=0.8, linewidths=2, zorder=5)
    
    ax.set_title(title, fontsize=11)
    ax.set_xlabel("X [mm]")
    ax.set_ylabel("Y [mm]")
    ax.set_aspect('equal', adjustable='box')
    
    return cf

# Create (100) and (110) contour maps for all depositions
for dep_name in DEPOSITIONS:
    if dep_name not in all_deposition_results:
        continue
    
    print(f"\nCreating (100) and (110) contour maps for {dep_name}")
    
    x_pos = all_deposition_results[dep_name]['x_coordinates']
    y_pos = all_deposition_results[dep_name]['y_coordinates']
    fwhm = all_deposition_results[dep_name]['fwhm']
    peakI = all_deposition_results[dep_name]['peak_intensity']
    peakA = all_deposition_results[dep_name]['peak_areas']
    N = x_pos.size
    
    fwhm = normalize_six(fwhm, N)
    peakI = normalize_six(peakI, N)
    peakA = normalize_six(peakA, N)
    
    planes_to_plot = [0, 1]  # (100) and (110)
    plane_names = [REGION_NAMES[i] for i in planes_to_plot]
    
    # ===== FWHM =====
    fig, axes = plt.subplots(1, 2, figsize=(12, 5), constrained_layout=True)
    for plot_idx, plane_idx in enumerate(planes_to_plot):
        plane_vals = fwhm[:, plane_idx]
        finite_vals = plane_vals[np.isfinite(plane_vals)]
        
        if finite_vals.size:
            pmin = float(finite_vals.min())
            pmax = float(finite_vals.max())
            if pmin == pmax:
                pmax = pmin + 1e-6
            plane_levels = np.linspace(pmin, pmax, 120)
        else:
            plane_levels = None
        
        cf = _contour_one_custom(axes[plot_idx], x_pos, y_pos, fwhm[:, plane_idx],
                                  title=f"{plane_names[plot_idx]}: FWHM (°)",
                                  cmap="coolwarm", levels=plane_levels)
        
        if cf is not None:
            fig.colorbar(cf, ax=axes[plot_idx], shrink=0.9, label="FWHM (°)")
    
    plt.suptitle(f"{dep_name}: FWHM Contour Map", y=1.00, fontsize=14)
    plt.show()
    
    # ===== Peak Intensity =====
    fig, axes = plt.subplots(1, 2, figsize=(12, 5), constrained_layout=True)
    for plot_idx, plane_idx in enumerate(planes_to_plot):
        plane_vals = peakI[:, plane_idx]
        finite_vals = plane_vals[np.isfinite(plane_vals)]
        
        if finite_vals.size:
            pmin = float(finite_vals.min())
            pmax = float(finite_vals.max())
            if pmin == pmax:
                pmax = pmin + 1e-6
            plane_levels = np.linspace(pmin, pmax, 120)
        else:
            plane_levels = None
        
        cf = _contour_one_custom(axes[plot_idx], x_pos, y_pos, peakI[:, plane_idx],
                                  title=f"{plane_names[plot_idx]}: Peak Intensity",
                                  cmap="coolwarm", levels=plane_levels)
        
        if cf is not None:
            fig.colorbar(cf, ax=axes[plot_idx], shrink=0.9, label="Peak Intensity [counts per second]")
    
    plt.suptitle(f"{dep_name}: Peak Intensity Contour Maps", y=1.00, fontsize=14)
    plt.show()
    
    # ===== Integrated Area =====
    fig, axes = plt.subplots(1, 2, figsize=(12, 5), constrained_layout=True)
    for plot_idx, plane_idx in enumerate(planes_to_plot):
        plane_vals = peakA[:, plane_idx]
        finite_vals = plane_vals[np.isfinite(plane_vals)]
        
        if finite_vals.size:
            pmin = float(finite_vals.min())
            pmax = float(finite_vals.max())
            if pmin == pmax:
                pmax = pmin + 1e-6
            plane_levels = np.linspace(pmin, pmax, 120)
        else:
            plane_levels = None
        
        cf = _contour_one_custom(axes[plot_idx], x_pos, y_pos, peakA[:, plane_idx],
                                  title=f"{plane_names[plot_idx]}: Integrated Area",
                                  cmap="coolwarm", levels=plane_levels)
        
        if cf is not None:
            fig.colorbar(cf, ax=axes[plot_idx], shrink=0.9, label="Integrated Area")
    
    plt.suptitle(f"{dep_name}: Integrated Area Contour Maps", y=1.00, fontsize=14)
    plt.show()
    
    # ===== Area Ratio =====
    area_100 = peakA[:, 0]
    area_110 = peakA[:, 1]
    
    with np.errstate(divide='ignore', invalid='ignore'):
        area_ratio = area_100 / area_110
    area_ratio[~np.isfinite(area_ratio)] = np.nan
    
    fig, ax = plt.subplots(1, 1, figsize=(8, 6), constrained_layout=True)
    
    finite_ratio = area_ratio[np.isfinite(area_ratio)]
    if finite_ratio.size > 0:
        rmin = float(finite_ratio.min())
        rmax = float(finite_ratio.max())
        if rmin == rmax:
            rmax = rmin + 1e-6
    
    cf = _contour_one_custom(ax, x_pos, y_pos, area_ratio,
                              title="Area Ratio: (100)/(110)",
                              cmap="coolwarm", levels=np.linspace(0.05, 0.25, 120))
    
    if cf is not None:
        fig.colorbar(cf, ax=ax, shrink=0.9, label="Area Ratio (100)/(110)")
    
    plt.suptitle(f"{dep_name}: Area Ratio Contour Map", y=1.00, fontsize=14)
    plt.show()

Creating (100) and (110) contour maps for 30-09-2025 Dep 1

Creating (100) and (110) contour maps for 22-09-2025 Dep 2

Creating (100) and (110) contour maps for 24-09-2025 Dep 3

Creating (100) and (110) contour maps for 07-10-2025 Dep 6

Creating (100) and (110) contour maps for 08-10-2025 Dep 7

Creating (100) and (110) contour maps for 09-10-2025 Dep 8

# Process each deposition and create distance plots
for dep_name in DEPOSITIONS:
    if dep_name not in all_deposition_results:
        print(f"\nSkipping {dep_name} - no results available")
        continue
    
    print(f"\n{'='*80}")
    print(f"Distance analysis for: {dep_name}")
    print(f"{'='*80}")
    
    # Get results for this deposition
    dep_data = all_deposition_results[dep_name]
    x_coordinates = dep_data['x_coordinates']
    y_coordinates = dep_data['y_coordinates']
    fwhm = dep_data['fwhm']
    peakI = dep_data['peak_intensity']
    peakA = dep_data['peak_areas']
    
    # Convert XRD coordinates to real-world coordinates
    # The x_coordinates from XRD need to have their sign flipped
    x_pos_real = -x_coordinates + wafer_center_x
    y_pos_real = y_coordinates + wafer_center_y
    z_pos_real = np.full_like(x_pos_real, wafer_center_z)
    
    # Calculate 3D distance from each point to target
    distances_3d = np.sqrt(
        (x_pos_real - target_x)**2 + 
        (y_pos_real - target_y)**2 + 
        (z_pos_real - target_z)**2
    )
    
    print(f"  3D Distance range: {distances_3d.min():.2f} - {distances_3d.max():.2f} mm")
    
    # Get data for (100) and (110) planes
    area_100 = peakA[:, 0]
    area_110 = peakA[:, 1]
    fwhm_100 = fwhm[:, 0]
    fwhm_110 = fwhm[:, 1]
    intensity_100 = peakI[:, 0]
    intensity_110 = peakI[:, 1]
    
    # Calculate area ratio
    with np.errstate(divide='ignore', invalid='ignore'):
        area_ratio = area_100 / area_110
    area_ratio[~np.isfinite(area_ratio)] = np.nan
    
    # ===============================================================
    # PLOT 1: FWHM vs 3D Distance
    # ===============================================================
    fig, axes = plt.subplots(1, 2, figsize=(14, 5), constrained_layout=True)
    
    # (100) FWHM
    ax = axes[0]
    mask_finite = np.isfinite(fwhm_100) & np.isfinite(distances_3d)
    ax.scatter(distances_3d[mask_finite], fwhm_100[mask_finite], 
               c='blue', s=50, alpha=0.7, edgecolors='black', linewidths=0.5)
    ax.set_xlabel('3D Distance to Target Center (mm)', fontsize=11)
    ax.set_ylabel('FWHM (°)', fontsize=11)
    ax.set_title('(100) Plane: FWHM vs 3D Distance', fontsize=12)
    ax.grid(True, alpha=0.3)
    
    if np.sum(mask_finite) > 3:
        z = np.polyfit(distances_3d[mask_finite], fwhm_100[mask_finite], 2)
        p = np.poly1d(z)
        x_trend = np.linspace(distances_3d[mask_finite].min(), distances_3d[mask_finite].max(), 100)
        ax.plot(x_trend, p(x_trend), 'r--', linewidth=2, alpha=0.7, label='Quadratic fit')
        ax.legend()
    
    # (110) FWHM
    ax = axes[1]
    mask_finite = np.isfinite(fwhm_110) & np.isfinite(distances_3d)
    ax.scatter(distances_3d[mask_finite], fwhm_110[mask_finite], 
               c='red', s=50, alpha=0.7, edgecolors='black', linewidths=0.5)
    ax.set_xlabel('3D Distance to Target Center (mm)', fontsize=11)
    ax.set_ylabel('FWHM (°)', fontsize=11)
    ax.set_title('(110) Plane: FWHM vs 3D Distance', fontsize=12)
    ax.grid(True, alpha=0.3)
    
    if np.sum(mask_finite) > 3:
        z = np.polyfit(distances_3d[mask_finite], fwhm_110[mask_finite], 2)
        p = np.poly1d(z)
        x_trend = np.linspace(distances_3d[mask_finite].min(), distances_3d[mask_finite].max(), 100)
        ax.plot(x_trend, p(x_trend), 'r--', linewidth=2, alpha=0.7, label='Quadratic fit')
        ax.legend()
    
    plt.suptitle(f'{dep_name}: FWHM by Distance to Target Center', y=1.00, fontsize=14)
    plt.show()
    
    # ===============================================================
    # PLOT 2: Peak Intensity vs 3D Distance
    # ===============================================================
    fig, axes = plt.subplots(1, 2, figsize=(14, 5), constrained_layout=True)
    
    # (100) Peak Intensity
    ax = axes[0]
    mask_finite = np.isfinite(intensity_100) & np.isfinite(distances_3d)
    ax.scatter(distances_3d[mask_finite], intensity_100[mask_finite], 
               c='blue', s=50, alpha=0.7, edgecolors='black', linewidths=0.5)
    ax.set_xlabel('3D Distance to Target Center (mm)', fontsize=11)
    ax.set_ylabel('Peak Intensity (counts per second)', fontsize=11)
    ax.set_title('(100) Plane: Peak Intensity vs 3D Distance', fontsize=12)
    ax.grid(True, alpha=0.3)
    
    if np.sum(mask_finite) > 3:
        z = np.polyfit(distances_3d[mask_finite], intensity_100[mask_finite], 2)
        p = np.poly1d(z)
        x_trend = np.linspace(distances_3d[mask_finite].min(), distances_3d[mask_finite].max(), 100)
        ax.plot(x_trend, p(x_trend), 'r--', linewidth=2, alpha=0.7, label='Quadratic fit')
        ax.legend()
    
    # (110) Peak Intensity
    ax = axes[1]
    mask_finite = np.isfinite(intensity_110) & np.isfinite(distances_3d)
    ax.scatter(distances_3d[mask_finite], intensity_110[mask_finite], 
               c='red', s=50, alpha=0.7, edgecolors='black', linewidths=0.5)
    ax.set_xlabel('3D Distance to Target Center (mm)', fontsize=11)
    ax.set_ylabel('Peak Intensity (counts per second)', fontsize=11)
    ax.set_title('(110) Plane: Peak Intensity vs 3D Distance', fontsize=12)
    ax.grid(True, alpha=0.3)
    
    if np.sum(mask_finite) > 3:
        z = np.polyfit(distances_3d[mask_finite], intensity_110[mask_finite], 2)
        p = np.poly1d(z)
        x_trend = np.linspace(distances_3d[mask_finite].min(), distances_3d[mask_finite].max(), 100)
        ax.plot(x_trend, p(x_trend), 'r--', linewidth=2, alpha=0.7, label='Quadratic fit')
        ax.legend()
    
    plt.suptitle(f'{dep_name}: Peak Intensity by Distance to Target Center', y=1.00, fontsize=14)
    plt.show()
    
    # ===============================================================
    # PLOT 3: Integrated Area vs 3D Distance
    # ===============================================================
    fig, axes = plt.subplots(1, 2, figsize=(14, 5), constrained_layout=True)
    
    # (100) Integrated Area
    ax = axes[0]
    mask_finite = np.isfinite(area_100) & np.isfinite(distances_3d)
    ax.scatter(distances_3d[mask_finite], area_100[mask_finite], 
               c='blue', s=50, alpha=0.7, edgecolors='black', linewidths=0.5)
    ax.set_xlabel('3D Distance to Target Center (mm)', fontsize=11)
    ax.set_ylabel('Integrated Area (a.u.)', fontsize=11)
    ax.set_title('(100) Plane: Integrated Area vs 3D Distance', fontsize=12)
    ax.grid(True, alpha=0.3)
    
    if np.sum(mask_finite) > 3:
        z = np.polyfit(distances_3d[mask_finite], area_100[mask_finite], 2)
        p = np.poly1d(z)
        x_trend = np.linspace(distances_3d[mask_finite].min(), distances_3d[mask_finite].max(), 100)
        ax.plot(x_trend, p(x_trend), 'r--', linewidth=2, alpha=0.7, label='Quadratic fit')
        ax.legend()
    
    # (110) Integrated Area
    ax = axes[1]
    mask_finite = np.isfinite(area_110) & np.isfinite(distances_3d)
    ax.scatter(distances_3d[mask_finite], area_110[mask_finite], 
               c='red', s=50, alpha=0.7, edgecolors='black', linewidths=0.5)
    ax.set_xlabel('3D Distance to Target Center (mm)', fontsize=11)
    ax.set_ylabel('Integrated Area (a.u.)', fontsize=11)
    ax.set_title('(110) Plane: Integrated Area vs 3D Distance', fontsize=12)
    ax.grid(True, alpha=0.3)
    
    if np.sum(mask_finite) > 3:
        z = np.polyfit(distances_3d[mask_finite], area_110[mask_finite], 2)
        p = np.poly1d(z)
        x_trend = np.linspace(distances_3d[mask_finite].min(), distances_3d[mask_finite].max(), 100)
        ax.plot(x_trend, p(x_trend), 'r--', linewidth=2, alpha=0.7, label='Quadratic fit')
        ax.legend()
    
    plt.suptitle(f'{dep_name}: Integrated Area by Distance to Target Center', y=1.00, fontsize=14)
    plt.show()
    
    # ===============================================================
    # PLOT 4: Area Ratio (100)/(110) vs 3D Distance
    # ===============================================================
    fig, ax = plt.subplots(1, 1, figsize=(8, 6), constrained_layout=True)
    
    mask_finite = np.isfinite(area_ratio) & np.isfinite(distances_3d)
    ax.scatter(distances_3d[mask_finite], area_ratio[mask_finite], 
               c='green', s=50, alpha=0.7, edgecolors='black', linewidths=0.5)
    ax.set_xlabel('3D Distance to Target Center (mm)', fontsize=11)
    ax.set_ylabel('Area Ratio (100)/(110)', fontsize=11)
    ax.set_title(f'{dep_name}: Area Ratio vs 3D Distance to Target', fontsize=12)
    ax.grid(True, alpha=0.3)
    
    if np.sum(mask_finite) > 3:
        z = np.polyfit(distances_3d[mask_finite], area_ratio[mask_finite], 2)
        p = np.poly1d(z)
        x_trend = np.linspace(distances_3d[mask_finite].min(), distances_3d[mask_finite].max(), 100)
        ax.plot(x_trend, p(x_trend), 'r--', linewidth=2, alpha=0.7, label='Quadratic fit')
        ax.legend()
    
    plt.show()
    
    # Print correlation statistics
    print(f"\nCorrelation coefficients (Pearson) with 3D distance:")
    for metric_name, data_100, data_110 in [
        ("FWHM", fwhm_100, fwhm_110),
        ("Peak Intensity", intensity_100, intensity_110),
        ("Integrated Area", area_100, area_110)
    ]:
        print(f"\n  {metric_name}:")
        for plane_name, data in [("(100)", data_100), ("(110)", data_110)]:
            mask = np.isfinite(data) & np.isfinite(distances_3d)
            if np.sum(mask) > 2:
                corr = np.corrcoef(distances_3d[mask], data[mask])[0, 1]
                print(f"    {plane_name}: r = {corr:.3f}")
    
    # Area ratio
    mask = np.isfinite(area_ratio) & np.isfinite(distances_3d)
    if np.sum(mask) > 2:
        corr = np.corrcoef(distances_3d[mask], area_ratio[mask])[0, 1]
        print(f"\n  Area Ratio (100)/(110): r = {corr:.3f}")

print(f"\n\n{'='*80}")
print("All depositions analyzed successfully!")
print(f"{'='*80}")

================================================================================
Distance analysis for: 30-09-2025 Dep 1
================================================================================
  3D Distance range: 158.64 - 224.99 mm

Correlation coefficients (Pearson) with 3D distance:

  FWHM:
    (100): r = 0.346
    (110): r = 0.340

  Peak Intensity:
    (100): r = 0.104
    (110): r = 0.120

  Integrated Area:
    (100): r = 0.062
    (110): r = 0.131

  Area Ratio (100)/(110): r = -0.041

================================================================================
Distance analysis for: 22-09-2025 Dep 2
================================================================================
  3D Distance range: 152.29 - 233.81 mm

Correlation coefficients (Pearson) with 3D distance:

  FWHM:
    (100): r = -0.830
    (110): r = -0.727

  Peak Intensity:
    (100): r = -0.849
    (110): r = -0.822

  Integrated Area:
    (100): r = -0.822
    (110): r = -0.843

  Area Ratio (100)/(110): r = -0.150

================================================================================
Distance analysis for: 24-09-2025 Dep 3
================================================================================
  3D Distance range: 158.64 - 224.99 mm

Correlation coefficients (Pearson) with 3D distance:

  FWHM:
    (100): r = -0.777
    (110): r = -0.494

  Peak Intensity:
    (100): r = -0.915
    (110): r = -0.939

  Integrated Area:
    (100): r = -0.879
    (110): r = -0.937

  Area Ratio (100)/(110): r = -0.601

================================================================================
Distance analysis for: 07-10-2025 Dep 6
================================================================================
  3D Distance range: 158.64 - 224.99 mm

Correlation coefficients (Pearson) with 3D distance:

  FWHM:
    (100): r = -0.737
    (110): r = -0.267

  Peak Intensity:
    (100): r = -0.889
    (110): r = -0.821

  Integrated Area:
    (100): r = -0.872
    (110): r = -0.796

  Area Ratio (100)/(110): r = 0.556

================================================================================
Distance analysis for: 08-10-2025 Dep 7
================================================================================
  3D Distance range: 158.64 - 224.99 mm

Correlation coefficients (Pearson) with 3D distance:

  FWHM:
    (100): r = -0.146
    (110): r = 0.099

  Peak Intensity:
    (100): r = -0.868
    (110): r = -0.786

  Integrated Area:
    (100): r = -0.843
    (110): r = -0.765

  Area Ratio (100)/(110): r = 0.448

================================================================================
Distance analysis for: 09-10-2025 Dep 8
================================================================================
  3D Distance range: 158.64 - 224.99 mm

Correlation coefficients (Pearson) with 3D distance:

  FWHM:
    (100): r = -0.750
    (110): r = -0.824

  Peak Intensity:
    (100): r = -0.904
    (110): r = -0.754

  Integrated Area:
    (100): r = -0.870
    (110): r = -0.782

  Area Ratio (100)/(110): r = 0.346


================================================================================
All depositions analyzed successfully!
================================================================================

XRD Analysis for All Depositions¶

Process All Depositions¶

Stacked Plots for All Depositions¶

Individual Plots with Smoothing and Peak Detection¶

Contour Maps - All FWHM¶

Contour Maps - All Peak Intensity¶

Contour Maps - (100) and (110) Planes: FWHM, Peak Intensity, Integrated Area and Ratio¶

Distance to Target Analysis for All Depositions¶